
pacman::p_load(tidyverse, sf, rnaturalearth, lme4, 
               gridExtra, raster, terra, geosphere, 
               nngeo, latex2exp, countrycode, 
               tmap, RColorBrewer, ggpubr, 
               spdep, spatialreg, sandwich,
               stargazer)



filter <- dplyr::filter
select <- dplyr::select

get_cluster_se <- function(model, cluster_var) {
  vcov_cl <- vcovCL(model, cluster = cluster_var)
  return(sqrt(diag(vcov_cl)))
}

#load data 
df<- read.csv("240617_global_method_scores_final.csv")

###########################################################################
#Table 2: Nomological validation
###########################################################################

lm1<- lm(data=df,
         vit_b_16 ~  log.rugged1k+river_binary)
lm2<- lm(data=df, 
         vit_b_16 ~ log.rugged1k+river_binary+
           wall_binary+log.distance+log.police.distance)
lm3 <- lm(data = df,  
          vit_b_16~ log.rugged1k+river_binary+
            wall_binary+log.distance+log.police.distance+bo_mean)
lm4<- lm(data=df, vit_b_16 ~ log.rugged1k+river_binary+
           log.distance+log.police.distance+wall_binary + bo_mean + avg_score_1.0km)
lm5<- lm(data=df, vit_b_16 ~ log.rugged1k+river_binary+
           log.distance+log.police.distance+wall_binary + bo_mean + avg_score_2.0km)
lm6<- lm(data=df, vit_b_16 ~ log.rugged1k+river_binary+
           log.distance+log.police.distance+wall_binary + bo_mean + avg_score_5.0km)
lm7<- lm(data=df, vit_b_16 ~ log.rugged1k+river_binary+
           log.distance+log.police.distance+wall_binary + bo_mean + avg_score_10.0km)

mods = list(lm1, lm2, lm3, lm4, lm5, lm6, lm7)
cluster_var <- ~ dyad
cluster_ses <- lapply(mods, get_cluster_se, cluster_var = cluster_var)

stargazer::stargazer(mods, 
                     se=cluster_ses, type='text',
                     omit.stat = c('ser', 'f', 'adj.rsq'))


###########################################################################
#Figures
################################################################################

legcorner<-st_as_sf(df, coords = c("xcoord", "ycoord"), crs = 4326)
legcorner <- st_make_valid(legcorner)

###########################################################################
#Figure 1
tmap_options(check.and.fix = TRUE)

raster_template <- raster(extent(legcorner), resolution = 0.5,
                          crs = 4326)
legrast <- rasterize(select(legcorner, vit_b_16), 
                     raster_template, 
                     field = 'vit_b_16', fun = mean)

world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf") %>%
  filter(geounit != "Antarctica") %>%
  filter(sovereignt != "Australia") %>%
  st_transform(crs = "+proj=eck6 +lon_0=0 +datum=WGS84")

# Plot with tmap
fig1<- tmap::tm_shape(world) +
  tm_fill("gray10") +
  tm_borders("black")+
  tm_shape(legrast) + 
  tm_raster(palette = rev(brewer.pal(8, "RdYlBu")),
            title = "Mean Border Legibility")+
  tm_layout(
    frame = FALSE,  # Remove frame
    inner.margins = c(0, 0, 0, 0),  # Remove inner margins
    bg.color = NA,  # Set background to transparent
    legend.position = c(0.7, 0.01),
    legend.title.size = 0.9,
    legend.text.size = 0.8
  )

#tmap_save(fig1, dpi=500, file='fig1.png')

###########################################################################
#Figure 2
world <- rnaturalearth::ne_countries(scale = "medium", returnclass = "sf")
world<- world %>% filter(geounit != "Antarctica") %>% st_set_crs(4326)
world <- st_make_valid(world)

ctry <- subset(world, sovereignt=="India" | sovereignt=="Pakistan")
states <- subset(ne_states(), admin=="India" | admin=="Pakistan")
kash <- subset(states, name=="Jammu and Kashmir" | name=="Azad Kashmir")

leg <- subset(legcorner, country1=="ind" & country2=="pak")
raster_template <- raster(extent(legcorner), resolution = 0.4,
                          crs = 4326)
points <- subset(legcorner, tile=="ind-pak-07089" | tile=="ind-pak-05073")

legrast_kash <- rasterize(dplyr::select(leg, vit_b_16), raster_template, 
                          field = 'vit_b_16', fun = mean)
bbox_new <- st_bbox(ctry) 
bbox_new[1] <- 60.84
bbox_new[2] <- 5.8
bbox_new[3] <- 90
bbox_new[4] <- 37.036670
bbox_new <- bbox_new %>%  
  st_as_sfc()

points$lab <- c("(A)","(B)")
points$cols <- c("white","black")
fig2<- tm_shape(ctry, bbox=bbox_new) +
  tm_fill("gray30") +
  tm_borders("black") + 
  tm_text("sovereignt") + 
  tm_shape(kash) +
  tm_fill("black") +
  tm_shape(legrast_kash) + tm_raster(palette=rev(brewer.pal(5,"RdYlBu")), 
                                     title="Mean Border Legibility") +
  tm_shape(points) + tm_dots(shape=7,size=0.2, col='cols' )+ 
  tm_text("lab",col="white", size=1,xmod=c(1.2,-1.2)) +
  tm_layout(legend.title.size = 1.2,
            legend.text.size = 1)

#tmap_save(fig2, dpi=500, file='fig2.png')

###########################################################################
#Figure 3

bbox_new <- st_bbox(world) 
bbox_new[1] <- 90
bbox_new[2] <- 0
bbox_new[3] <- 110
bbox_new[4] <- 30
bbox_new <- bbox_new %>%  
  st_as_sfc()

raster_template <- raster(extent(legcorner), 
                          resolution = 0.2,
                          crs = 4326)
legrast_sea <- rasterize(select(legcorner, vit_b_16), 
                         raster_template, 
                         field = 'vit_b_16', fun = mean)
temp<- subset(world, sovereignt!="Indonesia")

fig3<- tm_shape(temp, bbox=bbox_new) +
  tm_fill("gray20") +
  tm_borders("black") +
  tm_shape(legrast_sea) + 
  tm_raster(palette=rev(brewer.pal(8,"RdYlBu")),
            title="Mean Border Legibility") +
  tm_layout(legend.title.size = 1.2,
            legend.width = 2,
            legend.text.size = 1, 
            legend.title.color="black",
            legend.bg.color = TRUE)

#tmap_save(fig3, dpi=500, file='fig3.png')


###########################################################################
#SI

#Table S2: International borders with the highest legibility scores, on average

temp<- df %>% ungroup() %>%
  group_by(country1, country2) %>%
  summarise(vit_b_16_mean = mean(vit_b_16), 
            n = n()) %>%
  arrange(desc(vit_b_16_mean)) 

temp<- temp[1:10,]

tabs2<- xtable::xtable(temp, digits=3)


#Table S3: Nomological validation of border legibility (dyad random intercept model)

#With Dyad Random effect 

re1<- lmer(data=df, 
           vit_b_16 ~  log.rugged1k+river_binary+ (1|dyad))

re2<- lmer(data=df, 
           vit_b_16 ~ log.rugged1k+river_binary+wall_binary+log.distance+log.police.distance+
             (1|dyad))

re3<- lmer(data=df, 
           vit_b_16 ~ log.rugged1k+river_binary+wall_binary+log.distance+log.police.distance+bo_mean+
             (1|dyad))

re4<- lmer(data=df, 
           vit_b_16 ~ log.rugged1k+river_binary+wall_binary+log.distance+
             log.police.distance+bo_mean+avg_score_1.0km+
             (1|dyad))

re5<- lmer(data=df, 
           vit_b_16 ~ log.rugged1k+river_binary+wall_binary+log.distance+
             log.police.distance+bo_mean+avg_score_2.0km+
             (1|dyad))

re6<- lmer(data=df, 
           vit_b_16 ~ log.rugged1k+river_binary+wall_binary+log.distance+
             log.police.distance+bo_mean+avg_score_5.0km+
             (1|dyad))

re7<- lmer(data=df, 
           vit_b_16 ~ log.rugged1k+river_binary+wall_binary+log.distance+
             log.police.distance+bo_mean+avg_score_10.0km+
             (1|dyad))

mods = list(re1, re2, re3, re4, re5, re6, re7)

stargazer::stargazer(mods, 
                     omit.stat = c('f', 'ser', 'll', 'bic'),
                     omit = 'Constant')


# Table S4: Nomological validation of border legibility with lower resolution 

df<-st_as_sf(df, coords = c("xcoord", "ycoord"), crs = 4326)

raster_template <- raster(extent(df), resolution = 0.5,
                          crs = 4326)
legrast <- rasterize(select(df, vit_b_16), raster_template,
                     field = 'vit_b_16', fun = mean)
rugged_rast<- rasterize(select(df, rugged1k), raster_template,
                        field = 'rugged1k', fun = mean)
river_rast <- rasterize(select(df, river_binary), raster_template,
                        field = 'river_binary', fun = mean)
wall_rast <- rasterize(select(df, wall_binary), raster_template,
                       field = 'wall_binary', fun = mean)
pol_rast <- rasterize(select(df, log.police.distance), raster_template,
                      field = 'log.police.distance', fun = mean)
bc_rast <- rasterize(select(df, log.distance), raster_template,
                     field = 'log.distance', fun = mean)
bo_rast <- rasterize(select(df, bo_mean), raster_template,
                     field = 'bo_mean', fun = mean)


vit_b_16_values <- values(legrast)
ruggedness <- values(rugged_rast)
river_binary <- values(river_rast)  # Replace 'other_var_raster' with your raster
wall_binary <- values(wall_rast)
distance <- values(bc_rast)
poldistance <- values(pol_rast)
bo_mean<- values(bo_rast)

df2 <- data.frame(vit_b_16_values,ruggedness,
                  river_binary, 
                  wall_binary, 
                  distance, poldistance, bo_mean)


m1<- lm(df2, formula =vit_b_16_values ~ ruggedness+river_binary)
m2<- lm(df2, formula =vit_b_16_values ~ 
          ruggedness+river_binary+wall_binary+distance+poldistance)
m3<- lm(df2, formula =vit_b_16_values ~ 
          ruggedness+river_binary+wall_binary+distance+poldistance+bo_mean)


mods = list(m1,m2,m3)
ses = lapply(mods, function(x){sqrt(diag(vcovHC(x, "HC3")))})


stargazer::stargazer(mods, se = ses, 
                     omit.stat = c('ser', 'f'))